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Abstract. We study a generalized Dicke model, as recently realized in an atomic 
quantum gas experiment, describing the collective interaction of N two-level 
atoms with a single cavity mode. The model takes account of dissipation of 
the cavity field, and includes a non-linear atom-photon coupling, not present 
in the conventional Dicke model. We extend previous theoretical investigations 
of a semiclassical model by including all quantum effects and considering finite 
atom number N. Our results show good agreement between quantum expectation 
values and the semiclassical model as N is increased, but also show exotic 
behaviour for the corresponding quantum state as the non-linear atom-photon 
coupling is varied. 
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1. Introduction 

Fundamental models of many-body physics, originally developed as approximate 
theoretical descriptions of highly complex underlying systems, have in recent years 
been realized experimentally with an unprecedented degree of control thanks to 
advances in atomic physics and quantum optics. In particular, systems using 
ultracold atomic gases in optical lattices p] and systems based on cavity quantum 
electrodynamics (cavity QED) [5] have provided effective realizations (or "quantum 
simulations" [31 2]) of many-body models with a high degree of accuracy via exquisite 
experimental control and manipulation of atoms and electromagnetic fields. Such 
approaches also offer a flexible means of varying the effective model parameters, 
opening up exciting possibilities for exploring several frontiers of many-body physics 
and of strongly interacting quantum systems in general. 

Such realizations are particularly interesting in light of quantum phase transitions 
in many-body systems [5], as the flexibility of the schemes allows ready access to 
different phases and to the boundary (critical) regions between them. In particular, 
transitions between different phases can, for example, be observed in real time by the 
slow variation of one or more controllable system parameters (see, e.g., [2]). One can 
also consider more rapid "quench" dynamics by a sudden change of parameters, taking 
the system from one region of phase space to a distinctly different region over a short 
period of time (see, e.g., [B] and references therein) and thereby inducing complicated, 
non-equilibrium dynamics. In addition, signatures of the critical behaviour arc 
observable in the output channels from the systems; for example, in cavity QED 
systems, one typically finds dramatic changes in the properties of the output field of 
the cavity, such as the intensity and the photon counting statistics, and, for example, 
spontaneous symmetry breaking at a critical point may in fact be directly observable 
through homodyne detection of the field 6JJ 

Another interesting aspect of such effective realizations is that they often offer 
generalizations of the original model, sometimes in terms of non-trivial extensions of 
the Hamiltonian of the system, and also due to the intrinsically open (dissipative) 
nature of atomic and cavity QED systems, which in itself can lead to new phenomena. 
One particular example of this, which will be our focus in this work, is the effective 
realization of the Dicke model quantum phase transition using a superfluid gas 
in an optical cavity, recently demonstrated in [2]. Here, the two- level atoms of 
the original Dicke model are realized using resonant Raman transitions between 
a pair of discrete momentum states of a Bose-Einstein condensate (BEC). When 
comparing with the original Dicke model [7J [5], this realization offers a number 
of generalizations, such as the non-negligible decay of the optical cavity field, 
independently tunable interaction strengths for rotating and counter-rotating terms in 
the atom-field interaction Hamiltonian, and a new non-linear coupling term between 
the cavity photon number and the collective atomic inversion. The implications of 
these generalizations have been studied theoretically in great detail in j9[ [10] in the 
"thermodynamic limit" of a large number of atoms, where a semiclassical model can 
be employed. This analysis predicts an exceedingly rich dynamical phase diagram for 
the steady state of the system with fundamentally new behaviour and phases, such as 

f One should note that observing spontaneous symmetry breaking is typically extremely hard, as any 
external fluctuations or imperfections might destroy the necessary symmetry, but recent results have 
proved that atomic quantum gas experiments provide a promising route to observing such phenomena 
®. 
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a new supcrradiant phase, co-existence regions, and regimes with oscillatory long-time 
attr actors. 

The theoretical results of [9] HH] focus on parameters appropriate to the 
experiment in [2J. An analogous realization of the same model could also be based 
upon Raman transitions between a pair of electronic hyperfine ground states in alkali 
atoms, as first proposed and investigated theoretically in [11]. This leads to an identical 
effective Dicke model, but the scheme potentially offers access to a broader range of the 
effective parameters, and alternative ways of tuning them in order to access different 
regions of phase space. Further investigation of this scheme in [T2] , albeit focussing on 
the single-atom case, suggests that rubidium atoms coupled to a high finesse optical 
cavity should enable a good realization of the model. A many-atom cavity QED setup 
as realized recently in [13] is potentially of most immediate relevance to the current 
work. 

In particular, in this work we study this effective Dicke model including all 
quantum effects, and considering finite numbers of atoms. How this quantum model 
connects with the results for the thermodynamic limit considered previously is one 
aspect of our analysis, but, to further complement the work in [9j[T0], we also consider 
somewhat different parameter regimes, motivated by the possibility of realizing the 
model using the scheme discussed in the previous paragraph. We will see that one can 
achieve good acccess to the distinct regions of parameter space, with parameter ranges 
believed feasible based on the results from [12]. Furthermore, our findings indicate 
that one can observe phase transitions between a range of distinctly different phases 
through the variation of a single model parameter. 

2. Quantum model and semiclassical approximation 

The generalized Dicke model that we consider describes the interaction of N two-level 
atoms with a single mode of the electromagnetic field. The dynamics can be given in 
terms of a master equation for the density operator, p, of the atoms and cavity mode, 
as PI]: 



and I? is a super-operator defined through T>[0]p = OpO^ — l/20^0p — l/2pO^O, 
for any operator O. Parameters ujq and oj are the atomic and cavity frequencies, 
respectively, g is the linear interaction strength, and U is a non-linear coupling 
constant. Operators a and a* are the annihilation and creation operators for the cavity 
field, and {J 2 , J_, </+} are collective atomic operators satisfying angular momentum 
commutation relations: [J + , J_] = 2J Z , [J±,J Z ] = TJ±- Equation ([1]) reproduces the 
conventional equilibrium Dicke model if U = k = 0. 

This generalized Dicke model has been studied theoretically in [HI [TO] in the 
thermodynamic limit, N — > oo. In this limit one can derive a closed set of equations 
for the expectation values 



P=- 



i[H, p] + 2kT>[o\p, 



(1) 



where 




(2) 
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0= TV ' (3) 
TV 

Using equation (JTJ) and assuming factorization of operator products, i.e., 
(Jz,±{a + a*)) = (J z ,±) (a + a^) , / J_a^a) = (J_) (a^a) = (J_) | (a) | 2 , one can arrive 
at the following semiclassical equations of motion: 

a = — i (u) — in + Uj) a — ig(fi + /J*), 

/?= -i{Lo Q + U\a\ 2 )(3 + 2ig(a + a*) 1 , (4) 

7=i ff (a + o*)C8- J 8*). 

We note that both the quantum and semiclassical models conserve the total length of 
the spin. Physical initial conditions of course require that ( J x ) +(Jy) +( Jz) 2 < N 2 /A, 
where J± = J x ± iJ y . In [51 [TU] dynamical phase diagrams were mapped out by 
considering the long time attractors of equation ((4]). The semiclassical analysis predicts 
an exceedingly rich phase diagram, where the generalizations due to cavity decay 
and the non-linear coupling (U) leads to fundamentally new behaviour and phases, 
including a new superradiant phase, coexistence regions, and regimes of persistent 
oscillations for a and j3. Since one of our goals is to compare the semiclassical 
predictions with fully quantum results for finite N, we first summarize briefly some of 
the main results from [9j [10] . 

The steady states of equation (U) can be found analytically by requiring a = (3 = 
7 = 0. One finds that the semiclassical normal phase, {a = f3 = 0, 7 — —1/2}, and 
inverted phase {a = f3 = 0, 7 = +1/2} are always valid fixed points. In addition, two 
superradiant phases, referred to as superradiant A (SRA) and superradiant B (SRB), 
are possible. The SRA phase refers to solutions where (J y ) /N = i/2(/3 — j3*) = and 
the SRB phase to solutions where cj + U\a\ 2 = 0. We will assume w > throughout 
this paper, so that the latter condition is only possible for U < 0. The onset of 
superradiance from the normal/inverted phase was identified at critical values of the 
coupling constant g. The SRA phase comes into existence for g > g^ where + refers 
to the onset from the normal phase, — from the inverted phase, and 



The steady state values in this phase are 




g*{4u*-lP)-w UK 2 

U 2 (ujoU + 4g 2 ) ' [ba) 



U LO + (77 /1+ — IK 



y uj + U"f A —in 
Similarly, the SRB phase comes into existence if g > gB where 



9B = K u^-°(u/2ry u -°> (7) 
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and the steady state values are 



7 



B 



ir 



a 



B± 



= ±i 




(86) 



Since these results imply that multiple fixed points exists at any given point in 
parameter space, a further stability analysis of the different steady states is necessary 
to determine the phase diagram. However, one can immediately draw a few conclusions 
based on the above results. Physical solutions require 7 < 1/2, and of course 7 must 
be real. This means that for the two roots of the SRA solution j A± , the — solution 
can only exist for U < —2lj, and the + solution only for U < 2lu. Similarly, the SRB 
solution can only be valid if U < —2w (recall that we assume uiq > 0). The two critical 
points for U at U = ±2uj are found to mark the change in stability for the normal 
and inverted phases, the inverted phase being stable only if U < —2uj and the normal 
phase only if U < 2u>. We do not reproduce the scmiclassical stability analysis here, 
but refer to [10] for details. We note in particular that the region U > 2oj has no 
stable scmiclassical fixed points. It does yield a steady state solution for 7 given by 



as in the SRB phase, but the region has periodic orbits as long time attractors for a 
and f3. 

The results of the semiclassical steady state analysis can be summarized in a 
phase diagram for g and U. This is shown in figure [1] for parameters lu/(2it) = 1.0 
MHz, w /(2tt) = 0.05 MHz and k/(2tt) = 0.2 MHz (the specific values of which arc 
motivated by [12]). These results are analogous to the phase diagram in Figure 4 
of Ref. [9]. In particular, this diagram suggests that there is a possibility of taking 
the system through a whole range of different phases upon the variation of a single 
parameter, U. 

3. Numerical results for finite N 

In this section, we present numerical solutions of equation ([T]) as the parameter U is 
varied, and connect this to the semiclassical predictions discussed above. We focus 
our interest primarily on the steady state of equation (TTJ) , numerically approximated 
by an inverse power method iterated until an absolute tolerance of 10 -6 is achieved 
[14] . Alternatively, steady states could be approximated by integrating equation ([1]) 
for long times, but this approach is problematic due to the exceedingly long time- 
scales involved in certain regions of parameter space, as pointed out in [10] (see also 
section[3]4]bclow). The cavity mode Hilbcrt space was approximated with a truncated 
Fock space, with up to nine states for the largest value of g considered. 

We focus here on the parameters used in figure [T] i.e., oj/(2-k) = 1.0 MHz, 
w /(2tt) = 0.05 MHz, k/(2tt) = 0.2 MHz, U in the range -10.0 MHz < U/(2n) < 10.0 
MHz, and g/(2ir) up to 0.17 MHz. As for the number of atoms, we vary TV up to a 
maximum value of N = 90. 

A first observation that we would like to draw attention to, however, is that 
in regions of parameter space where the semiclassical analysis predicts co-existence 
phases, the size of the basin of attraction typically varies greatly for the individual 



U>2ui 



(9) 
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Figure 1. Schematic phase diagram based on the semiclassical approximation. 
The two superradiant phases are denoted SRA and SRB. JJ. refers to the normal 
phase {a, 0,j} = {0,0,-1/2}, f to the inverted phase {o,/3, 7} = {0,0, +1/2}, 
and "Limit Cycle" to the region U > 2lj, with periodic orbits as the semiclassical 
long time attractors for a and f3. The + sign indicates co-existence phases. The 
dashed lines show three horizontal cuts, at g/(2n) = {0.01,0.1,0.17} MHz, of the 
phase diagram that we will consider in the following sections. 

fixed points. If there are two, or more, fixed points with appreciably large basins 
of attraction compared to the characteristic size of quantum fluctuations, one would 
expect the steady state of the ensemble average dynamics given by equation (UJ to 
consist of a mixture of the solutions centered around the semiclassical fixed points. 
But, the semiclassical analysis does not predict the weighting of different fixed points, 
and we find this to in general depend on N and the other parameters. If a fixed 
point which is stable in the asymptotic limit, N — > 00, has a basin of attraction 
that is too small compared to quantum fluctuations for finite N, this fixed point 
will not contribute to the steady state of the quantum dynamics. We will often see 
this in our numerical results to follow. Another observation is that, for finite N, the 
normal (JJ-) and inverted (-ff) phases, |0,J| / it) = |0) cav ® \ J Z = z pN/2) sp i n , are in fact 
never the exact steady states of the quantum dynamics, as long as g > 0. Indeed, 
even in regions of parameter space where the semiclassical analysis predicts one of 
these phases as the unique stable fixed point, one finds a small, but non-zero, photon 
number in the cavity for finite N. This photon-production comes from the counter- 
rotating terms in equation (J5|), and might at first glance seem unphysical. We must 
here remember though, that equation ([1]) is an effective model that comes about via 
laser-pumped Raman transitions between distinct states of a multilevel system (such 
as an alkali atom or BEC), and the photon flux from the cavity can thus be explained 
by the energy input of the lasers involved in the effective realization. The normal and 
inverted phases are only approached in the limit N — > 00. 

With this in mind, we vary the parameter U along the horizontal cuts, shown 
as dashed lines in figure [TJ and compare the semiclassical predictions to steady state 
expectation values given by equation ([TJ . We also find it useful, as a means of shedding 
further light on the various phases and phase transitions, to look at how the reduced 
density operator for the cavity field, p ca ,v, and the Collective spin, yO S pin, change. This 
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Figure 2. Steady state expectation values, (Jz), (a^aj and (0), as U is varied 
along the lowest horizontal cut in figurc[l] {u)q,u>, K,g} = {0.05, 1.0, 0.2, 0.01} 2-7TX 
MHz: TV = 10 (dotted), TV = 50 (dashed), TV = 90 (solid). Semiclassical solutions 
are included as dashed green lines for comparison. 



can be visualized with, for example, the Wigner function and the atomic Q-function, 
respectively. Here the Wigner function for the cavity field is defined as 

W(a) = -tr \D^a)p cav D(a)(~lf a ] , (10) 

where a is a complex number and D(a) = exp(aa^ —a* a) is the coherent displacement 
operator. The atomic Q-function is defined by 

Q(v) = (v\p sP nM , (ii) 

where |ry) is the spin-coherent state, 

N //V\ 1/2 

W = (l + N 2 )- JV/2 E( n ) V n \n,N/2). (12) 



n=0 



Here, 77 = exp(i(j>) tan(0/2), with 9 and <f> spherical coordinates, and {\n,N/2)} are 
the symmetric Dicke states (see, e.g., [15]). 

3.1. Low g behaviour: g/(2n) = 0.01 MHz 

We start our discussion on the finite ./V results by considering the behaviour as U is 
varied at low g\ in particular, for g/(2n) = 0.01 MHz. This corresponds to the lowest 
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Figure 3. Contour plots of the Wigner function W(a), with a = (x + iy)/ y/2, for 
a series of values of U along the lowest horizontal cut in figure [l] {u>o,oj, K,g} = 
{0.05, 1.0, 0.2, 0.01} 2tt x MHz, N = 90. 
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Figure 4. Steady state atomic Q-function, Q(8,<f>), plotted on the Bloch 
sphere, for a series of values of U along the lowest horizontal cut in figure [T] 
{u ,u,k, g} = {0.05, 1.0,0.2, 0.01} 2tt x MHz, N = 90. The colour scheme goes 
from blue to red, indicating low to high concentration, respectively. Note that 
the individual subplots are not on the same color scale. 



dashed line in figure [T] There are no superradiant phases in this region of parameter 
space, but the semiclassical analysis suggests a transition from a co-existent JJ. + 
phase for U < — 2ui to a J} phase for — 2lj < U < 2lj. Furthermore, there is a transition 
to the "limit cycle" region for U > 2uj, where the semiclassical solutions for a(t) and 
/3(t) oscillate. It is not obvious from the semiclassical analysis how these transitions 
might manifest themselves in the quantum state of the system, i.e., the steady state 
of the quantum master equation ([1} . 

We first look, in figure [2l at the steady state expectation values (or "order 
parameters") (J 2 ) and (aJa) as U is varied along the cut g/(2n) = 0.01 MHz. We 
also show in the same figure the intensity correlation function of the cavity field, 
g( 2 \0) = (a^a^aa^ / (aJa) , in the steady state. We compare results for N = 10, 50 
and 90 atoms. The semiclassical solutions, as given in the previous section, are shown 
as green dashed lines wherever they apply, for comparison. As expected, the photon 
number in the cavity is small throughout, but there is a sharp increase at U = 2w. A 
pronounced increase is also seen in the collective atomic inversiop, (J z ), at this point, 
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with a shape that is characteristic of a second order phase transition. 

As for the region U < —2a;, the behaviour of the system is seen to vary greatly 
with N, and, for small N, also with U. Regarding the dependence on A, we point out 
that the energies of the states |0,ft) and |l,ft) are degenerate at U = — 2ojoA. The 
same is true for |0, ft) and |l,ft). Degeneracies in the Hamiltonian of the ground state 
and an excited state at a critical point is the trademark of equilibrium quantum phase 
transitions [5|. We are considering a dynamical steady state of a quantum master 
equation, but with a degeneracy in the spectrum at a value of U ~ N, certainly a 
significant dependence on A is not surprising. Note that this is not predicted in the 
thermodynamic limit, as the degeneracy is pushed out to U = — oo. One might wonder 
if there are any sharp changes to the steady state at this point of degeneracy in the 
spectrum. The expectation values in figure [2] do not indicate this. However, there is, 
in fact, a sharp transition in the quantum state of the spin as U crosses — 2o; A. This 
can be visualized using the atomic Q-function, as defined in the previous section and 
considered below. 

We study the reduced field and atomic states by looking at the cavity mode 
Wigner function and the atomic Q-function. Contour plots of the Wigner function for 
a series of values of U, with A = 90, are shown in figured] We see that the state of the 
field stays close to the vacuum for all values of U, with a symmetric Wigner function 
concentrated around the origin. The atomic Q-function is plotted on the Bloch sphere 
for the same values of U, and A = 90, in figure [4j Here, significant variation in the 
atomic state can be seen. Starting out as a ring shape for large, negative U, there 
is an increasing concentration around two points, approaching a two-peak structure 
as U is close to U = —9.0a;. As U crosses this point, the Q-function changes into a 
new two-peak structure, with a concentration around the north and south pole — this 
is what we would expect for a ft + ft co-existent phase. As U is increased further 
towards U = —2.0a;, the concentration shifts more and more towards the south pole, 
i.e., the ft phase, and at U = —2.5a;, no concentration around the inverted phase can 
be seen. The atomic state stays in close approximation to the ft phase throughout the 
interval —2a; < U < 2a;, but then there is a new sudden change around U = 2.0a;, 
where a ring like structure forms, and moves up towards the equator as U grows. 
This behaviour for the quantum state corresponds to the semiclassical predictions of 
a steady state value for (J z ), (J z ) /A = —uj/U, and periodic orbits for (J{ x ,y})- At 
large positive U, the state shares some similarities with the state at the corresponding 
negative U. 

3.2. s/(2tt) = 0.1 MHz 

Having studied the behaviour at low g, and seen some of the non-trivial finite- A effects 
that appear, we now move on to study the superradiant phases that can appear for 
larger g. In figure [5] we plot steady state expectation values as U is varied along 
g/(2n) = 0.1 MHz, i.e., the middle dashed line in figure [T] We compare three 
different atom numbers, A = 10, 30, and 50, and observe an increasing agreement 
with the semiclassical solutions as A grows. In particular, the results clearly point 
towards a first-order phase transition from a superradiant phase to the normal phase 
at U ~ —2.0a;. Comparing with the semiclassical predictions, we identify the phase for 
U < — 2oj as the SRB phase. Less clear is the second order phase transition from the 
normal to the SRA phase around U ~ 0.5a;. The semiclassical analysis also predicts 
a narrow normal phase region around 1.9a; < U < 2.0a;. For the finite A results the 
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Figure 5. Steady state expectation values, ( Jz), (at a} and g^ (0), as U is varied 
along the middle horizontal cut in figurc[T] {luq,uj, k, g} = {0.05, 1.0, 0.2, 0.1}27r x 
MHz: N = 10 (dotted), N = 30 (dashed), N = 50 (solid). Semiclassical solutions 
are included as dashed green lines for comparison. 



only evidence of this region is a slight dip in the photon number. This is perhaps not 
surprising given how narrow in parameter space this region is, and given the influence 
of quantum fluctuations. Then, finally, as U increases past U = 2.0a;, one enters the 
semiclassical limit cycle region, as in the low g case considered above. 

In figure[S]we show contour plots of the steady state Wigner function for the cavity 
mode, with atom number N = 50. We focus on values close to where the various phase 
transitions occur. At large negative U, we are in the superradiant phase, where the 
Wigner function has a two-peak structure. The two-peak structure becomes more 
clear as U is increased towards U = —2a;. Close to this value, the Wigner function 
rapidly changes into a single-peaked function centered at the origin, indicating the 
normal phase. The Wigner function starts splitting again as U crosses over beyond 
U = 0.5oj. The splitting is quite gradual, but one can see the Wigner function starting 
to deform around U = 0.4a;. At U = 1.5a;, the two-peak structure characteristic of 
a superradiant phase is clearly visible. At around U = 2.0a;, the Wigner function 
changes into an elongated shape around the origin, which we identify with the limit 
cycle region. In figure [7] we show the atomic Q-function at the same values of U. 
There is evidently a close correspondence between changes in the atomic state and 
changes in the state of the cavity field. 
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Figure 6. Contour plots of the Wigncr function W(a), with a = (x + iy)/ %/2, for 
a series of values of U along the middle horizontal cut in figure[T] {luo,lj, K,g} = 
{0.05, 1.0, 0.2, 0.1}2tt x MHz, N = 50. 
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Figure 7. Steady state atomic Q-function, Q(0,<j>), plotted on the Bloch 
sphere, for a series of values of U along the middle horizontal cut in figure [T] 
{£J ,w, K,g} = {0.05, 1.0, 0.2,0. 1}2tt x MHz, N = 50. The top two spheres 
are rotated 180° around the z-axis to better display their two-peak structure. 
The colour scheme goes from blue to red, indicating low to high concentration, 
respectively. The individual subplots are not on the same color scale. 



3.3. s/(27r) = 0.17 MHz 

Finally, we consider an even larger value for g. The top dashed line in figure [T] is 
at g/(2n) = 0.17 MHz. We show the variation in steady state expectation values, 
the Wigner function, and the atomic Q-function in figure |8l figure [9] and figure [TOl 
respectively. The results are quite similar to the results for g/(2n) = 0.1 MHz in the 
previous section, the main difference being that the transition at U = — 2.0lj is from 
the SRB phase to the SRA phase. In the Wigner function, this shows up as a rotation 
of the two peak structure in the x-y plane. The expectation values show reasonably 
good agreement with the semiclassical solutions even for N = 'SO. 



3.4. Steady state field correlation functions and power spectra 

An aspect, and indeed a concern, of some experimental significance is the large 
variation in the characteristic timescales of the evolution in different regions of 
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Figure 8. Steady state expectation values, (Jz), (aJa\ and g^ (0), as U is varied 
along the upper horizontal cut in figurcfTI {ujo, lj, k, g} = {0.05, 1.0, 0.2, 0.17} 2n X 
MHz: N = 10 (dotted), N = 30 (solid). Semiclassical solutions are included as 
dashed green lines for comparison. 



parameter space. The consequence of this for finite time experiments has been noted 
and discussed in some detail in [TO] • Here, we examine the characteristic timescales 
for the parameter ranges considered in the previous section by calculating the two-time 
amplitude correlation function for the cavity field, defined by 

C(t) := (at(i)a(0))-|(a) ss | 2 , 

where the steady state coherent intensity, | (a) ss | 2 , is subtracted for convenience. In 
figurcs fTTHroi we show the cavity field correlation function for the three different values 
of g previously considered, g/(2n) = {0.01,0.1,0.17} MHz, and for a series of values 
of U corresponding to those investigated earlier in the Wigner and Q-function plots. 

The correlation function plots illustrate clearly the wide range of dynamical time 
scales produced in the system through variation of the nonlinear coupling strength U, 
but we note in particular the exceedingly long-lived correlations for some parameter 
sets; in particular, for large negative U and the two larger values of g (i.e., U/uj = —10.0 
in figure [TO] and U/uj = —8.0 in figure ITO1) . The decay is not visible on the scale shown 
in the figures. These long-lived correlations can be related to the well-separated peaks 
in the atomic Q-functions shown for the corresponding parameters in figure [7] and 
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Figure 9. Contour plots of the Wigner function W(a), with a = (x + iy)/ %/2, for 
a series of values of U along the upper horizontal cut in figure [l] {luq,ui, «,<?} = 
{0.05, 1.0, 0.2, 0.17} 2tt x MHz, N = 30. 
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Figure 10. Steady state atomic Q-function, Q{9,<j>), plotted on the Bloch 
sphere, for a series of values of U along the upper horizontal cut in figure [T] 
{u ,u>,K,,g} = {0.05, 1.0, 0.2, 0.17} 2vr X MHz, N = 30. The first three spheres 
(starting from the top, traversing rows first) are rotated 180° around the z-axis. 
The colour scheme goes from blue to red, indicating low to high concentration, 
respectively. The individual subplots are not on the same color scale. 

figure [TU1 respectively, and, in particular, very weak coupling (or rate of "tunneling") 
between the two distinct states associated with these peaks (note that the cavity field 
is close to the vacuum state for these values of U). We also note long time scales 
involved close to the phase boundaries, U/ui ~ ±2.0, characteristic of critical slowing 
down. 

The Fourier transform of the correlation function, 

/oo 
dtC(t)exp{-iid), 
-oo 

gives the power spectrum of the cavity (output) field. We compute the spectrum 
numerically using a fast Fourier transform of the discrete time series shown in figures 
rrn ri3l This can be problematic due to the slow decay of the correlation function for 
some parameter sets, and we do not include spectra for those cases where the field- 
correlations have not decayed for times on the order of a few ms. In figure [14] and 

1 
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Figure 11. Field correlation function, C(t), for a series of U along the lower 
horizontal cut of figured] {^0,^,^,9} = {0.05, 1.0,0.2,0.01} 2tt x MHz, N = 90. 
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Figure 12. Field correlation function, C(t), for a series of U along the middle 
horizontal cut of figured] {^o,w, K,g} = {0.05, 1.0,0.2,0.1} 2vr X MHz, N = 50. 
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Figure 13. Field correlation function, C(i), for a series of U along the upper 
horizontal cut of figured] {^o,w, K,g} = {0.05, 1.0,0.2,0.17} 2tt x MHz, N = 30 



Generalized dissipative Dicke model 



20 



U/uj = -2.3 



to 




2000 
1000 


25 

20 

s 15 

to 10 

5 




U/uj = 0.4 



1 1 


, 1 


1 1 


u 


1 





U/uj = 2.0 





-0.08 0.00 0.08 
v jui 



-0.08 0.00 0.08 
v/uj 



Figure 14. Power spectra, S(ui), for a series of !7 along the middle horizontal 
cut of figurcQ] {tu ,u),K,g} = {0.05, 1.0, 0.2, 0.1} 2tt x MHz, N = 50. 



figure IT51 we show power spectra for the g/(2ir) = {0.1, 0.17} MHz cases, respectively. 
The power spectra provide an alternative (and directly measurable) means of studying 
the structure and dynamics of the system as U is varied; for example, very slow 
time scales in the dynamics are manifest as extremely sharp and pronounced spectral 
features. 
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Figure 15. Power spectra, S(ui), for a series of U along the upper horizontal cut 
of figure □ {uo,w,K,g} = {0.05, 1.0, 0.2, 0.17} 1% X MHz, TV = 30. 
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Figure 16. Logarithmic negativity in the steady state as U is varied along the 
horizontal cuts of figure [T] Top panel: {u)o,uj,K,g} = {0.05, 1.0, 0.2, 0.01} 2ir X 
MHz, N = 10 (dotted), N = 50 (dashed) and N = 90 (solid). Middle 
panel: {oj ,w,k, g} = {0.05, 1.0, 0.2, 0.1} 2tt X MHz, TV = 10 (dotted), N = 
30 (dash-dotted) and N = 50 (dashed). Bottom panel: {uo,u, K,g} = 
{0.05, 1.0, 0.2, 0.17} 2tt x MHz, N = 10 (dotted) and TV = 30 (dash-dotted). 



5.5. Field-spin entanglement 

Bipartite entanglement is believed to play a central role in quantum phase transitions, 
in the manner of divergence of correlation length at a critical point [T7| [T5] . This 
has been studied in the Dicke model in particular, both in terms of the bipartite 
entanglement between the cavity field and the collective spin, and in terms of 
entanglement between individual spins [THl QT] . 

Here we use the logarithmic negativity as a measure for the entanglement between 
the cavity mode and the collective atomic spin. This measure puts an upper bound on 
the amount of distillable entanglement between the two parts [21] . The logarithmic 
negativity is given by 

E N (p)=log 3 \\f?*\\i (13) 

where p r - A refers to the partial transpose of p over the subsystem A, and ||^||i = 
tr v X^X is the trace norm of X. In figure [TBI we show the logarithmic negativity in 
the steady state, as a function of U, for the three different cuts of the phase diagram 
(figured]) considered above: g/(2%) = {0.01, 0.1, 0.17} MHz. 
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Figure 17. Fidelity with the state 1 0, 0) , for the three different cuts of the 
phase diagram: {lu ,lu,k, g} = {0.05, 1.0, 0.2, 0.01} 2vr X MHz, TV = 90 (squares), 
{u> ,u>, K,g} = {0.05, 1.0,0.2,0.1} 2tt x MHz, N = 50 (crosses) and {w ,w, K, g} = 
{0.05, 1.0, 0.2,0.17} 2?r X MHz, N = 30 (diamonds). 



The results indicate a low degree of distillable entanglement for large N, except 
close to the supcrradiant phase transition around U ~ — 2lo (assuming g is large 
enough), in the SRA phase close to U = 2uj, and, interestingly, in the semiclassical 
limit cycle region, where U > 2.0w. For low g, the results add to what was observed 
in the expectation values considered earlier, namely a second order phase transition 
at U = 2lj, whereas the change in stability of the inverted phase at U — — 2u>, that 
was predicted by the semiclassical analysis, is not related to a phase transition for the 
quantum system. For higher g, however, the results point towards a hrst order phase 
transition at this value of U . 



4. Large \U\ behaviour 

The results so far point towards some similarities between large negative and positive 
U, especially for low g. The cavity photon number tends towards zero as \U\ grows, 
as does (J 2 ). Further investigation shows that for large \U\, and finite N, the atomic 
Q-function develops a distinctly ring-type structure, which moves towards the equator 
of the Bloch sphere with growing \U\, pointing towards a well-defined value for J z . 
Calculation of the purity of the steady state, tr[pg S ], further reveals the approach 
towards a pure state with growing \U\ and in turn towards a simplified dynamic for 
large \U\. In fact, investigation shows that the steady state approaches the state 
|0) cav ®\Jz = 0) in this limit, where \ J Z = j z ) is a J z eigenstate with eigenvalue j z . In 
figure[TT]we show the fidelity of the steady state with the state |0, 0) = |0) cav <8> | J z = 0) 
for 500 MHz < \U\/(2n) < 1000 MHz, where the fidelity of two quantum states p and 
g is defined as 



F(p,a) =tr 



(14) 
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5. Conclusions 



We have studied, by numerical means, a generalized Dicke model, where the long 
time attractor for the model is a dynamical quantum steady state, resulting from a 
balance between photon production, due to the counter-rotating terms of the Dicke 
model Hamiltonian, and decay of the cavity field. This generalized model has recently 
received attention due to its realization in a quantum gas experiment using an optical 
cavity. Other recent experimental demonstrations of self-organization for a large 
number of rubidium atoms coupled to a high-fincssc cavity also point towards a 
promising potential realization of the model. 

We have in particular focused on the consequences of a non-linear coupling 
between the cavity photon number and the collective atomic inversion, written 
as UJ z a^a/N, that is not present in the conventional Dicke model. If the non- 
linear coupling constant, U, becomes comparable to the cavity frequency, w, it can 
qualitatively change the steady state of the system. We have investigated this in 
detail, including quantum fluctuations and for a finite number of atoms, and seen how 
the steady state can be taken through a number of distinct phase transitions, which 
are observable, for example, as sharp changes in the cavity output field, as U is varied. 
The predictions of an approximate, semiclassical model are seen to be approached as 
the number of atoms is increased. 

For large values of the non-linear coupling, the atomic steady state approaches 
a state with a well-defined value for the inversion, J z , whereas J{ x ,y} become 
dclocalized — the latter corresponding to periodic orbits for the expectation values 
in the semiclassical approximation. In particular, as U tends to infinity, the steady 
state approaches the pure state |0) cav (& \J Z = 0) sp i n , i.e., the vacuum state for the 
cavity, while the atoms are in the maximally entangled J z = cigenstatc. For a finite 
number of atoms, N, this is also the case for U — > — oo, whereas in the thermodynamic 
limit, N — > oo, a mixture of J z = ±7V/2, eigenstates is predicted for large negative U. 
For finite N, this co-existence phase is observed for negative U < —2oj, but a sharp 
change in the atomic state appears at U ~ — 2wo-/V, corresponding to the development 
of a ring-type structure in the atomic Q'fimction, which eventually approaches the 
J z = eigenstate, as in the U oo case. 

The rich behaviour of the system, due to the presence of cavity decay and the 
non-linear atom-photon coupling, together with the readily observable cavity output 
field, suggests that various experimental implementations might offer observations of 
not only the conventional supcrradiant Dicke quantum phase transition, but also other 
novel quantum phases and phase-boundaries as well. We hope that this work might 
further stimulate experimental, as well as theoretical, investigations of the model. 
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